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Abstract. The numerical computation of the exponentiation of a real matrix 
has been intensively studied. The main objective of a good numerical method 
is to deal with round-off errors and computational cost. The situation is more 
complicated when dealing with interval matrices exponentiation: Indeed, the 
main problem will now be the dependency loss of the different occurrences 
of the variables due to interval evaluation, which may lead to so wide enclo- 
sures that they are useless. In this paper, the problem of computing a sharp 
enclosure of the interval matrix exponential is proved to be NP-hard. Then 
the scaling and squaring method is adapted to interval matrices and shown 
to drastically reduce the dependency loss w.r.t. the interval evaluation of the 
Taylor series. 



1. Introduction 

The exponentiation of a real matrix allows solving initial value problems (IVPs) 
for linear ordinary differential equations (ODEs): given A <E M Tlxn , the solution of 
the IVP defined by y'(t) = A y(t) and y(0) = yo is y(t) = exp(tA) yo, where for 
any M <E R nxn 

(1) exp( Af) : =£_ 

fc=0 

Linear ODE being met in many contexts, the numerical computation of the matrix 
exponential has been intensively studied (see |War77j . [BM89] , [MV03j . |Hig05| 
and references therein). While an approximate computation of leads to an 
approximate solution for the under liying IVP, interval analysis (see Section [2J offers 
a more rigorous framework: In most practical situations the parameters that define 
the linear ODE are known with some uncertainty. In this situation, one usually ends 
with an interval of matrices [A] = [A, A] := {A e R nx " : A < A < A} inside which 
the actual matrix A is known to be. Then, the rigorous enclosure of the solution 
will be obtained computing an interval matrix that encloses the exponentiation of 
the interval matrix [A]: 

(2) exp(L4]) := {exp(A) : A £ [A]}. 

The most obvious way of obtaining an interval enclosure of exp(A) is to evaluate 
the truncated Taylor series using interval arithmetic and to bound the remainder (cf. 
Subsection 14.11 for details). However, the next example shows that the truncated 
Taylor series is not well adapted to interval evaluation, even if no truncation of the 
series is performed. 

Example 1. Consider the interval of matrices A := [A, A] where 

(3) A:=(° Q \) ,A:=(° Q \) and Ait) := (j \ 
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Computing the formal expression of the exponential of the matrix A(t) for t £ 
[—3,-2], it can be proved that X_ < expQA]) < X with 



( 4 ) 

X 



_ , 1 i(l-e- 3 )\ ^ fl 0.316738 \ 
1 " e" 3 J ~ [p 0.0497871,) 
1 i(l-e- 2 )\ (1 0.432332\ 



2 e~' 2 I ~ \Q 0.135335/ 



where the lower and upper bounds cannot be improved, i. e. \X_, X] is the optimal 
enclosure o/exp([A]). 

Now, computing the interval Taylor series with an order of 10 (which is high 
enough to make the remainder insignificant) leads to T_< exp([A]) < T with 

1 -1.20912\ , = fl 1.95819\ 



(5) -~V° -6-25568J andT -y 6.4408 J" 

This is actually an enclosure of ([¥]), but a very pessimistic one. 

As shown by the previous example, even with high enough order for the expansion 
so the influence of the remainder is insignificant, the interval evaluation of the Taylor 
series computes very crude bounds on the exponential of an interval matrix. The 
reason of this bad behavior of the Taylor series interval evaluation is the dependency 
loss between the different occurrences of variable that occurs during the interval 
evaluation of an expression (cf. Section \2A\i . In general, one cannot expect to 
compute the optimal enclosure of ([2]): The NP-hardness of this problem is proved 
in Section [3] 

Two well known techniques can help decreasing the pessimism of the interval 
evaluation: First, centered forms can give rise to sharper enclosures than the natural 
interval evaluation for small enough interval inputs. Such a centered form for 
the matrix exponential was proposed in OM88a, OM88b, OM88c . However, this 
centered evaluation dedicated to the interval matrix exponentiation is quite complex 
and very difficult to follow or implement. Furthermore, there is an error in the proof 
of Proposition 10M of OM88aJ3 and some non justified assumptions in Section VII 
[OM88bj . 

The second technique consists of formally rewriting the expression so as to obtain 
a formula more suited to interval evaluation (usually decreasing the number of 
occurrences of variables) . For example, the evaluation of a polynomial in its Horner 
form is known to improve its interval evaluation [CG02] . It can be naturally applied 
to the Taylor series of the matrix exponential and was actually used in [OM88b 
to exponentiate the center matrix as required in the centered form. A proof of the 
correctness of the matrix exponential Taylor series Horner evaluation with rigorous 
bound on the truncation much simpler than the one given in [OM88bj is provided 
in Subsection |4?2l In Subsection l4.3[ we extend the well known scaling and squaring 
process (which consists of rewriting the Taylor series using the formula exp M = 
(exp M/2 S ) 2 ) to the exponentiation of interval of matrices. In addition of the 
usual benefits of this process, its use in conjunction with interval analysis allows an 
automatic control of the rounding errors. Furthermore, it is shown to drastically 



lr The proof of Proposition Proposition 10M of OM88a is claimed to be similar to the proof 
of Proposition 10 of |OM88a[ , However, the proof of Proposition 10 uses the fact that /([a;]) = 
Ui6[i]/W, which is valid only for scalar functions but not for vector- valued or matrix- valued 
functions, and thus cannot be extended to prove Proposition 10M which involves matrix-valued 
functions. 
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improve the dependency loss due to the interval evaluation, hence providing much 
more accurate and less expensive computations. As explained in Section[5]dedicated 
to experiments, the enclosure formula based on the scaling and squaring process 
is not only much simpler than the centered evaluation proposed in OM88b but it 
also provides sharper enclosures. 

2. Interval Analysis 

Interval analysis (IA) is a modern branch of numerical analysis that was born in 
the 60's. It consists of computing with intervals of reals instead of reals, providing a 
framework for handling uncertainties and verified computations (see [M0066, AH74, 
INeu901 IJKDWOI] and [Kea96j for a survey). 

2.1. Intervals, interval vectors and interval matrices. An interval is a con- 
nected subset of R. Intervals are denoted by bracketed symbols, e.g. [x] CI, When 
no confusion is possible, lower an upper bounds of an interval [x] are denoted by 
x £ K and with x < x, i.e. [x] = \x, x] — {x £ M : x < x < x}. Furthermore, 
a real number x will be identified with the degenerated interval [x,x\. 

There are two equivalent ways of defining interval matrices. On the one hand, 
being given two matrices A < A £ R nxm (where the inequality is defined compo- 
nentwise), an interval of matrices is obtained by considering 

(6) [A] := {A e R nxrn : A < A < A}. 

On the other hand, being given intervals [ffly], a matrix of intervals is obtained by 
considering 

(7) [A] := {A G R nxm : V^ £ {1, . ..n},Vj 6 {1, . ..m} 7 a tJ G [a tJ ]}. 

These two definitions are obviously equivalent following the notational convention 
A = (cty), A = (aij) and [ay] = [0^,0^], and will be used indifferently. 

Interval vectors are defined similarly to interval matrices as either intervals of 
vectors or vectors of intervals. 

The magnitude of an interval [x] is |[x]| := max{|a:|, \x\}. The magnitude of 
an interval matrix is the real matrix formed of the magnitude of the entries of 
the interval matrix, i.e. \[A]\ := (|[tty]|)ii- ^ ne infinite norm will be considered 
in the rest of the paper. The norm of an interval matrix is the maximum of the 
norms of the real matrices included in this interval matrix. It is easily computed 
as I p]|| = || P]| ||. It obviously satisfies [A] C [B] implies ||[A]|| < ||[B]||. 

The interval hull of a set of real is the smallest interval that contains this set. It 
is denoted by □. For example □{(), 1} = [0, 1]. The interval hull is defined similarly 
for sets of vectors and sets of matrices. 

2.2. Interval Arithmetic. Operations o G {+, x, — , -=-} are extended to intervals 
in the following way: 

(8) \x,x] o [y,y] ■= {xoy : x G [x,x],y G [y,y]}- 

The division is defined for intervals [y, y] that do not contain zero. Note that 
unary elementary functions like exp, In, sin, etc., can also be extended to intervals 
similarly. All these elementary interval extension form the interval arithmetic (IA). 
As real numbers are identified to degenerated intervals, the IA actually generalizes 
the real arithmetic, and mixed operations like 1 + [1, 2] = [2, 3] are interpreted using 
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The I A lacks some important properties verified by its real counterpart: It is 
not a field anymore as interval addition and interval multiplication has no inverse 
in general, while distributivity is not valid anymore (instead a subdistributivity 
law holds in the form + [z]) C [a;] [y] + [x] [z]). On the other hand, interval 

operations are inclusion-increasing, i.e. [x] C [x'] and [y] C [y'\ imply [x] o [y] C 
[x']o[y']. 

Rounded Computations. As real numbers are approximately represented by floating 
point numbers |Gol91j , the IA cannot match the definition (j5J) exactly. In order to 
preserve the inclusion property, the IA has to be implemented using an outward 
rounding. For example, [l,3]/[2, 2] = [0.5,1.5] while both 0.5 and 1.5 cannot be 
exactly represented with floating point numbers. Therefore, the computed result 
will be [0.5+ 1.5 + ] where 0.5~ (respectively 1.5 + ) is a floating point number smaller 
than 0.5 (respectively bigger than 1.5). Of course, a good implementation will 
return the greatest floating point number smaller than 0.5 and the smallest floating 
point number greater than 1.5. Among other implementations of IA, we can cite the 
C/C++ li braries PROFIL/BIAS [K nu94j an d Gaol |Gou06j . the Matlab toolbox 
INTLAB }Rum99] and Mathematica |Wol08j . 

2.3. Interval Evaluation of an Expression. The natural usage of the IA is 
to evaluate an expression for interval arguments. The fundamental theorem of 
interval analysis (cf. |Moo66| ) allows explaining the interpretation of this interval 
evaluation. Its proof is classical but is reproduced here. 

Theorem 1. Let E and F be either K or M" or R nxn and [x] G IE. Consider a 
real function f : [a;] — ► F and an interval function [/] : I[x] — > IF, where I[x] is 
the set of all intervals included in [x], i.e. l[x] = {[y] G IE : [y] C [x]}. Suppose 
furthermore that both 

(1) For all xe[x], f{x) G [f](x) 

(2) [/] is inclusion-increasing inl[x]. 

Then, [f]([x])D{f(x):xe[x]}. 

Proof. For all x G [x] we have f(x) G [f](x) by (1), and because [x,x] C [x], (2) 
implies [f](x) C[f]([x]).D " □ 

Let us illustrate the usage of the fundamental theorem of IA on a simple example, 
which can be trivially generalized to arbitrary expressions. Consider the expression 
x + xy. When evaluated for degenerated interval arguments, it gives rise to [x, x] + 
[x, x] [y,y] — [x + xy, x + xy] 3 x + xy. Furthermore, it is inclusion-increasing as it is 
compound of inclusion-increasing operations. Therefore, the fundamental theorem 
of I A proves that [x] + [x] [y] D {x + xy : x G [x] , y G [y]}. As another example, the 
interval evaluation of the expression of the interval matrix/matrix product [A] [B] 

(9) (MPOtf = E Mi^] 

k 

gives rise to the inclusion [A] [B] D {AB : A G [A],B G [B]}. 

When the expression evaluated for interval arguments contains only one oc- 
currence of each variable, the computed enclosure of the range is optimal. As a 
consequence, the expression (j9|) is optimal since only one occurrence of each vari- 
able is involved in the expression of each entry. Thus the more accurate statement 
[A] [B] = D{AB : A G [A], B G [B]} actually holds. Note that [A] [B] + 1 {AB : A G 
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[A] , B G [B] } since the product [^4] [B] actually contains matrices that are not the 
product of matrices from [^4] and [B], but [A] [B] is the smallest interval matrix 
that contains {AB : A G [A],B e [B]}. 

However, the interval evaluation of expression that contains several occurrences 
of some variable is not optimal anymore in general. In this case, some overesti- 
mation generally occurs which can dramatically decrease the usefulness of interval 
evaluation. 

2.4. Over-estimation of Interval Evaluation. When an expression contains sev- 
eral occurrences of some variables its interval evaluation generally gives rise to 
a pessimistic enclosure of the range. For example, the evaluation of x + xy for 
the arguments [x] — [0,1] and [y] = [—1,0] gives rise to the enclosure [—1,1] of 
{x + xy : x £ [x],y G [y]} while the evaluation of a;(l + y) for the same interval 
arguments gives rise to the better enclosure [0, 1] of the same range (the latter en- 
closure being optimal since the expression x(l + y) contains only one occurrence of 
each variables). This overestimation is the consequence of the loss of correlation be- 
tween different occurrences of the same variables when the expression is evaluated 
for interval arguments. 

While [A] [B] = D{AB : A G [A],B G [B]}, the interval evaluation of [A] [A], 
which encloses {A 2 : A G [A]}, is not optimal in general since several occurrences of 
some entries of [A] appear in each expression of the entries of [A] [A] . An algorithm 
for the computation of D{A 2 : A G [A]} which can be evaluated with a number 
of interval operations that is polynomial w.r.t. the dimension of [A] was proposed 
in |KKMN0"5] . However, it was proved that no such polynomial algorithm exists 
for the computation of Dj^ 3 : A G [A]} unless P=NP, i.e. the computation of 
□ {^ 3 : A G [A]} is NP-hard (cf. [KKMN05Q . 

The situation is even worth than this: Even computing an enclosure of {A 3 : 
A G [A]} for a fixed precision is NP-hard. The notion of e-accuracy of an enclosure 
is introduced to formalize this problem (see e.g. (Gag'85, KLR.K98]). The following 
definition is adapted to sets of matrices. 

Definition 1. Let A C R" xm be a set of matrices, [A] = DA G IM. nxm , and 
consider an interval enclosure [B] of A (which obviously satisfies [B] "D [A]). The 
interval enclosure [B] is said e-accurate if 

(10) max! max [a- — b ti \ , max \a,ij — 6«| } < e 

ij ij 

Thus, an e-accuracy enclosure of a set of matrices is e-accurate for each entry. 
Although it is not stated in [KKMN05], the proof presented there also shows that 
the computation of an e-accuracy enclosure of {A 3 : A G [A]} is NP-hard. 

3. Computational Complexity of the e-AccuRATE Interval Matrix 

Exponentiation 

Computing e-accurate interval enclosures of the range of a multivariate polyno- 
mial / : R" — ► M is NP-hard (cf. |Gag85| and Theorem 3.1 in [KLRK98Q . Even if 
one restricts its attention to bilinear functions, the computation of e-accurate en- 
closures of their range remains NP-hard (cf. Theorem 5.5 in [KLRK98]). Note that 
if one fixes the dimension of the problems, then the computation of these e-accurate 
enclosures is not NP-hard anymore, hence showing that the NP-hardness is linked 
to the growth of the problem dimension. It is not a surprise that computing an 
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e-accurate enclosure of the interval matrix exponential is NP-hard, but this result 
remains to be proved. 

Theorem 2. For every e > 0, computing an e-accurate enclosure o/exp(L4]) is 
NP-hard. 

Proof. We prove that the e-accurate enclosure of the range of a bilinear function, 
which is NP-hard by Theorem 5.5 in [KLRK98 , reduces to the e-accurate enclosure 



of the interval matrix exponential. 



(11) 



A := 



( 


T 

X 





°\ 








B 














y 


V o 








o J 



Let B 6 



and A := 



and [x],[y] € IR". Define 



( 


x T 





°\ 








B 














y 


V o 








o J 



which are obviously computed in polynomial time from B, [x] and [y]. We now 
prove that an e-accurate enclosure of the exponentiation of [A] gives rise to an e- 
accurate enclosure of the range of the image of [x] and [y] by the function x T B y, 
which will conclude the proof. Let 
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be such that A € [A], that is equivalently x € [x] and ye [y] . One can check easily 
that A is nilpotent: 
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Thus (exp(A)) 1 2n+2 = \ x T B y. As a consequence, 

(14) {(exp(A)) 1 2n+2 : A e [A]} = {x T B y : x e [x],y e [x]} 

and the entry (l,2n + 2) of an e-accurate enclosure of exp([A]) is an e-accurate 
enclosure of the image of [x] and [y] by the function x T B y. □ □ 

4. Polynomial Time Algorithms for the Enclosure of the Interval 

Matrix Exponential 

This section presents three expressions dedicated to the enclosure of the expo- 
nential of an interval matrix: The naive interval evaluation of the Taylor series, the 
interval evaluation of the Taylor series following the Horner scheme and the interval 
evaluation of the series following the scaling ans squaring process. 

4.1. Taylor Series. The naive interval evaluation of the truncated Taylor series 
for interval matrices exponential is now presented including a rigorous bound on 
the truncation error. The bound used here is the same as in in |OM88bj . Let us 
define for K + 2 > |p]|| 

[f]([A},K):= I + [A] + \[Af + ... + ^[A] K 
{ ) [T]{[A],K):= [f]([A],K) + [K]([A},K), 
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where the interval remainder [1Z] ([A], K) is 
(16) [R]([A],K) :=p(\\[A}\\,K) [-E,E] with p(a,K) 



a K+1 



and E E M. nxn has all its entries equal to 1. We provide now a new proof that 
[T] ( [A] , K) is an enclosure of {exp(A) : A E [A]} which is much simpler than the 
one provided in |OM88bj and which will be used in the rest of the paper. The 
following lemmas will allow applying the fundamental theorem of interval analysis 
to expressions that include [1Z] ( . , K) . 

Lemma 1. For a fixed positive integer K , the interval matrix operator [1Z]( . , K) 
is inclusion-increasing inside {[A] E IR nxrl : < K + 2}. 

Proof. Let [A], [B] E M nx ™ such that < K + 2 and [A] C [B]. Then, \ \[A}\\ < 

\\[B}\ \ which implies | [A] 1 1 < K + 2. Furthermore as p(a, K) is obviously increasing 
with respect to a, we have p(||[A]||, K) < p(\\[B]\\, K). Finally, as [~E,E] is 
centered on the null matrix, we have 

(17) P(\\[A]\\,K) [-E,E}Cp(\\[B]\\,K) [-E, E], 

which concludes the proof. □ □ 

The next lemma a direct consequence the well known upper bound on the trun- 
cation error of the exponential series. 

Lemma 2. Let A e R nx " and K G N such that K + 2 > \\A\\. Then exp(A) <E 
[T](A,K). 



Proof. Suppose that exp(A) ^ [T](^4, if), i.e. there exist i,j E {1, ...n} such that 
(18) (ex P (A)) 2j ,i (E^r) r + 



k=0 



This obviously implies 



(19) \(ex P (A)).-C£^) i] \> p(\\A\lK). 

Therefore ||exp(A) — X)^o4rll > Pdl^ll*-^) holds, which contradicts the well 
known bound on the truncation error for the exponential series (see e.g. [OM88b ). 
Eventually exp(^) E [T](A,K) has to hold. □ □ 

Theorem[3]below states that [T]([A], K) is an enclosure of exp([A]). It was stated 
in [OM88bj but proved with different arguments in |OM88aj . [OM88bj . Note that 
the usage of the fundamental theorem of interval analysis allows us to provide a 
proof much simpler than the one proposed in [OM88a , |OM88bj . 

Theorem 3. Let [A] E lR nxn and K E N such that K + 2 > \\[A}\\. Then 
exp([A])C[T]([A],K). 

Proof. First, by Lemma [1] [1Z]( . , K ) is inclusion-increasing, and therefore so is 
[T]( . , K ) because it is compounded of inclusion-increasing operators. Second, by 
Lemma[2](V/l E [A]) exp(A) E [T](A,K). Therefore, one can apply the fundamen- 
tal theorem of interval analysis to conclude the proof. □ □ 
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Example 2. Consider the interval of matrices [A] defined in Example^ Theorem 
[3 with K — 16 gives rise to the following enclosure of exp(L4]): 

'l + [-9 x 1(T 7 , 9 x 1(T 7 ] [-1.2092, 1.9582]\ 
[-9 x 10~ 7 , 9 x 1(T 7 ] [-6.2557, 6.4409]/ ' 

Higher order for the expansion do not improve the entries (1,2) and (2,2) anymore. 



(20) 



4.2. Horner scheme. The Horner evaluation of a real polynomial improves both 
the computation cost and the stability (see e.g. |Knu97j ). When an interval evalua- 
tion is computed, the Horner evaluation furthermore improves the effect of the loss 
of correlation (see |CG02| ). It is therefore natural to evaluate (fTS"]) using a Horner 
scheme: 



(21) 



[H]([A},K):= I+[A][l+ [ 4[l+ [ 4[ ... (l + l§ 
[H]([A],K):= [H}([A},K) + [R]([A},K). 



(22) 



Lemma 3. Let A e R nxn and K € N such that K + 2 > |L4||. Then exp(A) e 
[H](A,K). 

Proof. As interval operations are evaluated with real arguments, the Horner scheme 
can be expanded exactly leading to [H](A,K) = [T](A, K). As a consequence, 
[H]{A, K) = [T](A, K) and Lemma M concludes the proof. □ □ 

Theorem 4. Let [A] € IR" X " and K € N such that K + 2 > \\[A}\\. Then 
exp([A])e[H]([A],K). 

Proof. As a consequence of LemmaQ] [H] ( . , K) is compounded of inclusion- increasing 
operator and is hence inclusion increasing while Lemma [3] shows that exp(^4) G 
[W](^4, K). Therefore, one can use the fundamental theorem of interval analysis to 
conclude the proof. □ □ 

Example 3. Consider the interval of matrices [A] defined in Example]]^ Theorem 
2] with K = 16 then to the following enclosure o/exp([/l]): 

'1 + [-1.1 x 10~ 6 , 1.1 x 10~ 6 ] [-0.0706, 0.7352]\ 
[-1.1 x 10~ 6 , 1.1 x 10~ 6 ] [-1.2056, 1.2117y ' 

This enclosure is sharper than the one computed using the Taylor series: as it was 
forseen, the Horner evaluation actually improves the loss of dependency introduced 
by the interval evaluation in the expression of the Taylor expansion of the matrix 
exponential. 

4.3. Scaling and squaring process. The scaling and squaring process is one of 
the most efficient way to compute a real matrix exponential. It consists of first 
computing exp(^4/2 L ) and then squaring L times the resulting matrix: 

(23) exp(A) = (cxp(^/2 L )) 2i . 

Therefore, one first has to compute exp(^4/2 L ). This computation is actually much 
easier than exp(A) because ||A/2 L || can be made much smaller than 1. Usually, 
Pade approximations are used to compute exp(^4/2 L ). However, this technique has 
not been extended to interval matrices, hence we propose here to use the Horner 
evaluation of the Taylor series instead. Therefore, we propose the following operator 
for the enclosure of an interval matrix exponential: Let K and L be such that 
(K + 2)2 L >\\[A]\\ and define 
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(24) [S]([A],L,K):=([H]([A]/2 L ,K)) . 

The exponentiation in (|24|) is of course computed with L successive interval matrix 
square operations. 

Theorem 5. Let [A] G IE nx " and K, L G N such that (K + 2) 2 L > \\[A]\\. Then 
exp([A])C[S]([A},L,K). 

Proof. By Theorem^ we have exp(A/2 L ) G [H] (L4]/2 L , K) for an arbitrary ^4 G 
[A]. The interval evaluation [M] 21 * of an arbitrary interval matrix [M] encloses 
{NP L : M G [M]}, and therefore [S]([A],L,K) B exp{A/2 L f L . This concludes the 
proof as this holds for an arbitrary A G [A]. □ □ 

Example 4. Consider the interval of matrices [A] defined in Example^ Theorem 
[5| with L = 10 and K = 10 /eads to i/je following enclosure o/exp([A]): 

, , /l + [-5.7 x 10- 13 , 9.1 x 10- 13 ] [0.3165, 0.4325]\ 

^ > \ [-2.4 x 10~ 19 ,2.4 x 10- 19 ] [0.0496, 0.1355],) ' 

This enclosure is much sharper than the two previously computed using the Taylor 
series (cf. Example^) and its Horner evaluation (cf. Example 0). It is also 
very close to the optimal enclosure (j4]). The computation cost for L and K is 
approximately the same as the Horner scheme with order L + K . 



5. Experiments 

In this section, we compare the interval Horner evaluation of the truncated Taylor 
series versus the interval scaling and squaring method. The direct Taylor series is 
not presented as it is similar but always worth that its Horner evaluation. In order 
to compare these two enclosures, we use the width of these interval enclosure: Let 
wid [M] be the real matrix formed of the widths of the entries of [M]. We will use 
the 1 1 wid [M]\\ as a quality measure of the enclosure [M]. Experimentations have 
been carried out using Mathematica [W0IO8J . Subsection 15.11 presents a detailed 
study of the exponentiation of one real matrix, while Subsection 15.21 deal with an 
interval matrix exponentiation. 

As explained in introduction, the comparisons presented in this section do not 
include the interval matrix enclosure method proposed in OM88a, OM88b. OM88cJ. 
However, this method is based on a interval Horner evaluation and thus cannot give 
rise to better enclosures than the interval Horner evaluation of the center matrix, 
which is of poor quality as demonstrated above. 

5.1. Interval exponentiation of a real matrix. In this section, we consider the 
matrix A defined by 

/ 131 19 18\ 
(26) A := -390 56 54 

\-387 57 52/ 

proposed in [BM89]. This matrix is difficult to exponentiate because it has signifi- 
cant eigenvalue separation and a poorly conditioned eigenvector set. 
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||wid [H](A,K)|| 
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FIGURE 1. Plots of ||wid [H] ( A, K) \ | w.r.t. K, for two different 
precisions: Plain line for p = 110 digits precision, and dashed line 
for p — 120 digits precision. 



The matrix A is too difficult to exponentiate using the Horner evaluation of the 
truncated interval Taylor series: As \\A\\ = 500 the Taylor series requires an expan- 
sion of order greater than 502. Computing [H](A,K) using double precision does 
not provide any meaningful enclosure. Figure [T] shows the quality of the enclosure 
obtained for different orders ranging from 1500 to 1800 and two different precisions 
for computations (the Mathematica [W0IO8] arbitrary precision interval arithmetic 
was used) . It shows that no meaningful enclosure is obtained for precision less than 
p = 110 digits or order less than K — 1550 using the Horner interval evaluation of 
the Taylor expansion. 

The interval scaling and squaring formula gives rise to ||wid [S] (A, 12, 12) 1 1 « 
7.2 x 10 -6 computed using the standard double precision arithmetic. This can be 
improved using a decomposition A = PMP^ 1 where M is easier to exponentiate. 
Then, one can compute expA = P exp(P _1 AP)P (where P _1 has to be rig- 
orously enclosed in order to maintain the rigorousness of the process). Using the 
Shur-dccomposition, we obtain 



5.2. Interval exponentiation of an interval matrix. In order to compare the 
different methods, we will use 0.1^4 which is simpler to exponentiate. We have expo- 
nentiated [A e ] := 0AA+[— e, e] for various values of e inside [10~ 16 , 1] and the results 
are plotted on Figure [2j The three plain gray curves represent 1 1 wid [i?]([/l e ], K)\\ 
for K = 150, K = 160 and K = 170 (increasing K improves the enclosure, until 
K = 170 above which no significantly improvement is shown). The black curve 
represents ||wid [,5]([A e ], 10, 10)||. The dashed line represents 



(27) 



\\wid P[S]{P 



1 AP, 12,12) P- 1 ]] w 7.2 x 10 
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(28) 



I wid D{exp A^ , exp A e }\\, 
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FIGURE 2. Log- log graphics for the comparison of the quality for 
different enclosing methods. 



which is a lower bound of ||widDexp(L4 e ])||. Each plot show three phases: The first 
phase shows flat plots, then ||wid ( • )|| increases linearljQ w.r.t. e before eventually 
exponentially increasing. 

During the flat phase, the rounding errors represent the main contribution to 
the final width of the enclosure. Thus, decreasing e does not decrease the width of 
the final enclosure. 

The linear phase is the most interesting. For these values of e, the width of 
□ exp([A e ]) growth linearly because the contribution of quadratic terms are negli- 
gible. On the other hand, the interval evaluations [i?]([A e ], K) and [/S]([A e ], 10, 10) 
are pessimistic, but it is well known that the pessimism of interval evaluation grows 
linearly w.r.t. the width of the interval arguments. Thus, the computed enclosure 
show a linear growth w.r.t. e which are approximately 

(29) |wid [H]([A f ],K)\\ « 1.17 x 10~ 4 + 2.86 x 10 10 e 

(30) ||wid[S](p4 e ], 10, 10)|| « 1.80 x 10~ 9 + 8.59 x 10 3 e. 

This cleary shows how smaller is the pessimism introduced by the interval scaling 
and squaring process. 

Finally, both the interval Horner evaluation and the interval scaling and squaring 
process show an exponential growth when e is too large. The lower bound repre- 
sented by the dashed line also shows a exponential growth, which proves that this 
is inherent to the exponentiation of an interval matrix. For such e, some matri- 
ces inside [A e ] eventually see some of their eigenvalues becoming positive, leading 
to some exponential divergence of the underlying dynamical system, which is also 
observed on the matrices exponential. 



The linear plot displayed within the log-log scale indicates a polynomial behavior, the poly- 
nomial degree being fixed by the slope in the log-log representation. Here, the slope is one inside 
the log-log plot and thus so is the degree of the polynomial. 
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